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ABSTRACT 

This paper examines a number of analytic models which can describe the 
gravitationally stabilized fusion reactor of the Sun. An analytical multipa- 
rameter model is shown to reproduce various structure parameters such as 
density, mass, pressure and temperature throughout the solar core, which 
have been obtained by solving numerically the system of solar structure 
differential equations. For simplicity, numerical and analytical results are 
discussed in detail for a two-parameter model. 

1 INTRODUCTION 

The numerical approach to the problem of studying solar structure is to go 
for numerical solutions of the underlying system of differential equations (A. 
Noels, R. Papy, and F. Remy 1993). Even for a simple main-sequence star 
in hydrostatic equilibrium, like the Sun, at least four nonlinear differential 
equations are to be dealt with to obtain a detailed picture of the run of 
physical variables throughout the star. Another approach is to search for 
analytical solutions of these differential equations under justified physical 
simplifications (Haubold and Mathai 1992). 

In this paper we start with numerical data given in Sears (1964) for various 
analytical models which will fit the data on density in the gravitationally 



stabilized solar core. Let r be the distance from the center of the Sun to an 
arbitrary point in its interior, R the solar radius, p(r) the density at r, p c 
the density at the center and y = r/R. For the range < y < 0.3 consider 
the following models for the density p(r). 



u = ^ = 1 - 4.94y + 6.67y 2 - 2.73y 3 (1.1) 

Pc 

u = 1 - Ay + 2y 2 + 2y 3 - y 4 (1.2) 

u = (i-V^)a-y 3 r (i-3) 

u = (l-y 3/2 ) 16 (1.4) 

u = (l-v^)(l-l/ 3 ) 64 (l-2/) (1-5) 

u = (1-y 1 - 48 ) 14 (1.6) 

u = (l-y L48 ) 13 (1.7) 

u = (1-y 1 - 28 ) 10 (1.8) 

The following table gives y, the numerically obtained data on u = p(r)/p c 



from Sear's (1964) and the estimated values under models (1.1) to (1.8). 
Table 1 density distribution 



V 

y=R 


(Sears 


, 1964) Model 


Model 


Model 


Model 






(1.1) 


(1.2) 


(1.3) 


(1.4) 


0.0864 


0.6519 


0.6213 


0.6720 


0.6795 


0.6626 


0.1153 


0.5253 


0.5149 


0.5690 


0.5987 


0.5283 


0.1441 


0.3856 


0.4185 


0.4710 


0.5722 


0.4065 


0.1873 


0.2810 


0.2908 


0.3340 


0.3720 


0.2588 


0.2161 


0.1994 


0.2164 


0.2490 


0.2796 


0.1837 


0.2450 


0.1424 


0.1499 


0.1700 


0.1957 


0.1264 


0.2882 


0.0962 


0.0649 


0.0560 


0.0985 


0.0679 



Table 1 continued 



Model 


Model 


Model 


Model 


(1.5) 


(1.6) 


(1.7) 


(1.8) 


0.6208 


0.6849 


0.7037 


0.6418 


0.5297 


0.5573 


0.5811 


0.5229 


0.4435 


0.4404 


0.4669 


0.4179 


0.3023 


0.2928 


0.3196 


0.2885 


0.2192 


0.2159 


0.2409 


0.2202 


0.1478 


0.1551 


0.1772 


0.1649 


0.0701 


0.0885 


0.1053 


0.1033 



Model (1.1) results from a least square fit of a third degree polynomial to 
the data of Sears (1964). Model (1.2) comes from successive eliminations. 
Even though these models can adequately describe the data in the range 
< y < 0.3 the equation u=0 has real roots in (0,1) besides the root at 
y — 1. Hence these cannot be used as models for the density in [0,1]. From 
models (1.3) to (1.8) it is evident that a model of the type 

u = (l-y a ^(l-y a *) b *...(l-y a *) b *, (1.9) 

where > 0, fej = 0, 1, 2, . . . , i — 1, . . . , k can be an excellent fit to the data. 
For simplicity we will study a model of the type 

u = (l-y s )\ (1.10) 

where 5 > and 7 a positive integer. Illustration will be given for the case 
5 = 1.28 and 7 = 10. 

2 A Solar Model 

We assume that rotation and magnetic fields have no impact on the internal 
structure of the Sun and the gravitational force directed inward and the gas 
pressure force directed outward keep the gas sphere in hydrostatic equilib- 
rium. The ratio of gas pressure and radiation pressure increases towards the 
center of the Sun but never exceeds 10~ 3 . Thus radiation pressure can be 



neglected for the purpose of considering hydrostatic equilibrium. Assume 
that the density p{r) varies from p c at r = to zero at r = R in the following 
fashion: 

f _ 

p(r) = p c [l — (— ) ] 7 , 5 > 0, 7 a positive integer. (2.1) 
R 

Then the distribution of mass in terms of radius is given by 

= 47rr 2 p(r), 



dM{r) 2 



so that 



M(r) = 4tt dtt 2 p(t) 
Jo 

= Anp c J^dtt 2 [l-(^) 
= ^p c jr(-ir(l)R* r H dzz 



r/R 

m5+2 



3 K R J u <T <S ' v iT ; ' v ; 

where 2-^1 is a Gauss' hypergeometric function, see for example, Mathai 
(1993). M(r) gives the total mass at r = R, that is, 

4tc Pc R 3 „ . 3 3 
M(i2) = 2 F 1 (- 7 ,-;- + l; 1) 



4-Kp c R 3 7! 



(2.3) 



3 (f + l)(f + 2)...(f + 7)' 

which is evaluated by using the formula 

, , . T(c)T(c-a-b) . . 

2 f 1 ( n ,6;c ;1) = r ; e ^ n)r(e _ 6 j . (2.4) 

From (2.2) and (2.3) we can get a formula for estimating p c also. That is, 

_ 3M(fl)(f + l)(f + 2)...(f + 7 ) 
Pc " "4^ 7! ' { } 

The numerical values for the Sun from Sears (1964) are 

p c = 158gcm~ 3 
M(R) = 1.991 x 10 33 # (2.6) 
R = 6.96 x 10 10 cm. 



Substituting these in (2.5) we have 



(i+m+f--(i^) = 112m (2 . 7) 

7! 

Hence a pair of (5, 7) which satisfy (2.7) is a good choice. We have many- 
choices for (5, 7) even if 7 is kept as a positive integer. Keeping 7 a positive 
integer, for simplicity, it is found that (5 = 1.28,7 = 10) is a convenient 
choice in the sense that these values give good estimates for M(R) and R 
also, close to the realistic values and the estimates for the density are not far 
off. 

Note that (2.7) is obtained by using the readings in (2.6). From (2.7) one 
can get a set of positive values for 5 and 7. These are given here for some 
integer values of 7: 
Table 2 7 and 5 



7 


2 


3 


4 


5 


6 


7 


5 


0.2225 


0.4412 


0.6265 


0.7802 


0.9098 


1.0211 




8 


9 


10 


11 


12 


13 




1.1182 


1.2043 


1.2814 


1.3512 


1.4149 


1.4735 




15 


16 


17 


18 


19 


20 




1.5780 


1.6251 


1.6692 


1.7108 


1.75 


1.7872 



14 



From (2.2) and (2.3) the relative mass is given by 

M(r) _ (f + l)(f + 2)...(f + 7) r 3 3.3 r s 

M{R) ~ ~\ { R ] 2Fl( ~ 7 ' 6' 5 + 1; { R ] } - (2 ' 8) 

This is computed for 5 = 1.28, 7 = 10 and compared with Sear's (1964) 
numerical results in the following table 3. 



Table 3 mass distribution 



y=R 


M(r) 
M(R) 

(analytic model) 


M(r) 
M(R) 

(Sears 1964) 


0.0864 


0.0533 


0.05 


0.1153 


0.1106 


0.1 


0.1441 


0.1866 


0.2 


0.1873 


0.3249 


0.3 


0.2161 


0.4241 


0.4 


0.2450 


0.5225 


0.5 


0.2882 


0.6576 


0.6 



Assuming that the pressure P(r) at the center of the Sun is P c and that 
at the surface is zero we have 

P(r)=Pc _ G r dt mm, (2 . 9) 

Jo t z 

where G is the gravitational constant. That is, 

P{r) = P c _^ p ^]T<-^/ dM r[l-(4)'] 7 

5 ^ m\ Jo (l + m) 1 V J 



2 2 1 ,r 

- + m-- + m + l-(- 



x 



2i r i(-7,T + m-- + m + l-(-) s ). (2.10) 



Assuming that the pressure at the surface is zero, that is P(R) = 0, we 
obtain 

AnG i (- 7 ) ro 1 

f c - —^-PcK — " X 



Then 



5 2 lc ml (| + m )(| + m ) 

7j 

+ m + 1)(| + m + 2) . . . (| + m + 7) ' 



(2.11) 



2 F 1 (- 1 ,- + m;- + m + l;(^) s ) 



(2.12) 



The above expansion can also be written as a double sum. In this case 

l)m {-l)n 



P(r) = P t 



AitO 7 7 

. e e 



6 



x 



m=0 n=0 
5 ) m( 5) m( 5) m+n 



ml 



ni 



y s + l)rra(x + l)m(i + l)m+n 



"(-)'" 






n 






K R J . 





(2.13) 



This double sum can be written in terms of a Kampe de Feriet's function. 
That is, 



P{r) = Pc~ ^GpyFlZZ 
and then assuming that P(R)=0, 



p c = 1>kGpIp?fI^ a; 



(*) 4 



3 2 



2._<v 3 2. 
i • S ' S ■ ' 



where the Kampe de Feriet's series is defined by the following: 

pp:q:k \(l 
r:m:n 



r x\ |(%>MM : ( c fc) 
\V) K a r):(/3 m ):(7„) 



(2.14) 



(2.15) 



(2.16) 



[n^ =1 (%) m / +n / 




nj =1 (6j) m ' 




nf = i(cj)„'_ 




[n^ =1 (aj) m / +n / 








n™ =1 (7,)n' 


m'ln'l 



E E 



For a discussion of Kampe de Feriet's function see Srivastava and Karlsson 
(1985). Note that (2.14) is a polynomial since 7 is a positive integer and 
hence convergence conditions do not arise in (2.14). The form in (2.10) is 
the most appropriate for computational purposes and the computations given 
later are done by using (2.10). 

The temperature is given by the equation of state of the perfect gas 



(2.17) 



where // is the mean molecular weight, k is Boltzmann's constant and Na is 
Avogadro's number. Let 



g(r) 



1 £ (-7), 

m=0 ""■ \8 



5 2 ir'n m! (§ + m)(f+m) 



x 



7! 



(f + 771+ !)...(§ +771 + 7) 



\ra<5+2 



R' 



s Fi(-7, - + m; - + m + 1; 



Then 



P(r) = 47rGp^ 2 ^(r) 
and the temperature is then given by 



(2.18) 
(2.19) 



T(r) 



kN, 



-AnGp c R 2 



[1 " 



(2.20) 



In the following table we tabulate g(r) for 5 = 1.28 and 7 = 10 from which 
we can compute the pressure for the solar core of the sum ranging from 



< £ < 0.3. 



Table 4 pressure = const. xg(r), temperature = const. xg(r)/u, 



u 



1 - (j^Y , for 5 = 1.28 and 7 = 10 



V = Ti 

0.0864 
0.1153 
0.1441 
0.1873 
0.2161 
0.2450 
0.2882 



g(r) 

1.7911 x 10~ 3 
1.4035 x 10~ 3 
1.0546 x 10~ 3 
6.4117 x 10- 4 
4.4168 x 10- 4 
2.9482 x 10- 4 
1.5279 x 10~ 4 



0.6418 
0.5228 
0.4179 
0.2884 
0.2202 
0.1649 
0.1033 



g( r ) 

u 

2.7907 x 10~ 3 
2.6842 x 10~ 3 
2.5234 x 10~ 3 
2.2227 x 10~ 3 
2.0055 x 10~ 3 
1.7876 x 10~ 3 
1.4788 x 10~ 3 



Hence the proportional decrease in temperature is: 
Table 5 temperature distribution 



oeais yLvyjHij 


g(r) 




9(r) 


[P(r)/Pc] 1/2 




\p(r)/ P c] 


0.8789 


0.9428 




0.8987 


0.9275 


0.9316 




0.8819 


8828 


0.9219 




8385 


0.8938 


0.9612 




0.8956 


0.8911 


0.9709 




0.8992 


0.9000 


0.7644 




0.8656 


Table 6 pressure 


= const. xg(r), for 8 


= 1.28 


T 

y = R 


g(r) 






0.0864 


3.1681 x 10~ 


-3 




0.1153 


2.7142 x 10" 


-3 




0.1441 


2.2664 x 10- 


-3 




0.1873 


1.7286 x 10- 


-3 




0.2161 


1.4425 x 10- 


-3 




0.2450 


1.2012 x 10- 


-3 




0.2882 


0.9129 x 10- 


-3 





3 ENERGY GENERATION 

From Mathai and Haubold (1988) we have the net release of thermonuclear 
energy per gram per second given by 

e[p(r),T(r)} = eo [p(r)]"[T(r)] m 



kN A p{r) 



for the parameter n and m determined by the proton-proton chain and the 
CNO cycle. Since the g(r) of table 5 is proportional to P(r) we have 

e[p(r),T(r)}=kMr)} n - m [9(r)] m , (3.2) 



where k\ is a constant. The equation of energy conservation is 

47rr 2 p(r)e[p(r),T(r)], (3.3) 



dr 

where L(r) denotes the luminosity of the Sun. Hence the total luminosity is 
given by 

rR 

L = 4tt / drr 2 p(r)e\p(r),T(r)] (3.4) 
J o 

= Ane ^Y j\rr'[p(r)r n - m [P(r)r (3.5) 
for the model in (3.2). Substituting for P(r) from (2.10) one has 

L = (47T6 ) m+1 (^) m pr n+1 (^) m X 



L 



R drAl - (^F (1+n - m) [*(7,5) - 0(^)] m , (3.6) 



where 



mi =o "-i- u ' ""i;va 

2 i r i(-7,T T + l + m i; (^)) (3.7) 



and 

W7)= E 



# m ^o m i ! b + m i)(j + m i 

-l)mi 1 7 



X 



, (ll =o ™i ! (f + m 1 )(| + m 1 )(| + m 1 + l)---(| + m 1 +7)' 

(3.8) 

Note that | | < 1. If m and n are not assumed to be positive integers 
then the general binomial expansion can be used but one has to check the 
convergence conditions in the final sum. Assuming that m and n are positive 
integers we get polynomials from the two factors in the integral in (3.6). For 
a general positive integer m we get higher powers of $(-u) but $(w) itself is a 
double sum. Hence the expression becomes complicated. For m=l and n a 
positive integer one can get some simpler representations. Hence we consider 
this case. Then we have 

L = (4 7 re ) 2 (^)pr 2 (^) J R 3 [<f (5,7)/o - hi (3.9) 



where 



and 



i = f 1 duu 2 [i-u s r = (n7) 

Jo 



r(f) 



5 r(| + n 7 +l) 



(3.10) 



h 



duu 2 [l-u 5 ] n ^(j){u) 

■l)mi (-7)t 



X2 7 7 



, g)mi(g)mi 



mi=0m2 =o m i ! m s! (| + l) mi (f + l) mi 



X 



(!) 



gjm 1 +m2 



{g + l)mi+m 2 



d UU 2 [l - ^]n 7 U mi5W+2 . (3.11) 



But 



Jo 



(n 7 )! r(f + mi+m 2 ) 
5 r(| + n 7 + l + m 1 + m 2 ) 



(n 7 )! 



r(f) 



5 r(§ + n 7 + l) 



x 



( <5 )mi+m2 



(t + n 7 + l) mi+ma 



(3.12) 



Substituting (3.12) in (3.11) and writing the resulting expansions as Kampe 
de Ferier's function we have 



5 r(*) 
J ' = 6 <"'' )! r(| + n -, + 1) 



x 



F 2:3:l 
r 2:2:0 



(1) 



2 5. 2 3. 
|+l,f+n 7 +l:|+l,|+l: 



The luminosity at any given point r is available from (3.3). That is, 

,GR 2 



L(r)= (47re ) 



m+l 



\m ^,m+Tt+l ( ft \m 



(3.13) 



(3.14) 



dtt 2 



t 1 - <l' s 


7(l+n— m) - 







1) 



We will evaluate this for m=l. Then the two integrals to be evaluated are 



h = 



dtt 1 



r' 



n-y 



and 



That is, 



f 

Jo 



dtf 



R' 



™7 f 



h = R 3 J T/R duu 2 [l -w* 5 ]" 7 
= R s g (-n7)« rr/R 

a=0 



duu 2+aS 







= ^3 g (-^7). (r/i?) 3+a5 

a=0 



3 + aS 



r 3 "7 

-E 

° a=0 



(-^7)c 



2^1 



a! (f + l). 
3 3 



r V 



(3.15) 



jj>3 ^ ^ ^ ( l)mi ( 7)t»2 ( n l)m 3 ^ 
mi=0m 2 =0m 3 =0 m l' 



mi+m 2 



(i) 



m 2 ! m 3 ! 

5+(mi+m,2+m3)5 



(3.16) 



(| + l) 5 [§ + m 1 + m 2 + m 3 



5i? 2 



mi +m2 

7 7 



'-7)mi (-7) 



m 2 l < 5Jmi+m 2 



X 



. 5 Jmi+rre 2 



r V 



1711+TO2 



X 



(5 + l)mi+m 2 

2 Fi ^-717, ^ + m 1 + m 2 ; ^ + 1 + mi + m 2 ; {^j j . 

Note that (3.16) can also be written in terms of a Kampe de Feriet's function 
of three variables. In that case the variables are 1, (r/R) s , (r/R) s , respec- 
tively. Note also that (3.15) and (3.16) are polynomials. Hence for m=l, n 
a positive integer 



L(r) = (4tt 



,2 / \ „n+ 2 



S 2 



Pc 



\kN A 



[*(7,5)/2-/ 3 ], 



(3.17) 



where I2 and I3 are given in (3.15) and (3.16) respectively. A particular case 
of (3.17) is also available from Haubold and Mathai (1992). 
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